## One dimensional IRT model
## Ordered response options

## N survey respondents
## J policy topics asked
## K treatment/endorsement options (including control) 

data {
 dim.vill <- dim(villcov)
 dim.dist <- dim(distcov)
 dim.indiv <- dim(indivcov)
}

model{
  for(i in 1:N){ ## loop over respondents
  
    x[i] ~ dnorm(delta.cov10 %*% indivcov[i, ] + delta.village0[village[i]], .02)   
  
    for (j in 1:J){ ## loop over questions            

      s[i,j] ~ dnorm(lambda.cov10[endorser[i,j], ] %*% indivcov[i, ] + lambda.village0[village[i],endorser[i,j]], pow(omega[endorser[i,j]], -2))
      
      p[i,j,1] <- phi(alpha[j,1] - beta0[j]*(x[i] + s[i,j]))
      p[i,j,2] <- phi(alpha[j,2] - beta0[j]*(x[i] + s[i,j])) - (p[i,j,1])
      p[i,j,3] <- phi(alpha[j,3] - beta0[j]*(x[i] + s[i,j])) - (p[i,j,1] + p[i,j,2])
      p[i,j,4] <- phi(alpha[j,4] - beta0[j]*(x[i] + s[i,j])) - (p[i,j,1] + p[i,j,2] + p[i,j,3])
      p[i,j,5] <- 1 - (p[i,j,1] + p[i,j,2] + p[i,j,3] + p[i,j,4])
      
      Y[i,j] ~ dcat(p[i,j,1:5])
      
      ## In-sample average prediction rates
      correct[i,j] <- p[i,j,Y[i,j]]
    }}
  
  for (j in 1:J){ beta0[j] ~ dlnorm(0, 1) } ## log-normal prior for discrimination parameter  
  beta <- beta0 * sd.x ## If sd.x is large, that forces beta0 to be small
  
  for (j in 1:J){ for (l in 1:4){ alpha0[j,l] ~ dnorm(0, .02) }}
  for (j in 1:J){ alpha[j,1:4] <- sort(alpha0[j,]) } ## Cutpoints

  for (l in 1:dim.indiv[2]) {
    delta.cov10[l] ~ dnorm(0, .01)
  }
  delta.cov1 <- delta.cov10 / sd.x

  for (l in 1:dim.indiv[2]) {
    lambda.cov10[1, l] <- 0
    for (k in 2:K){
      lambda.cov10[k, l] ~ dnorm(0, .01)
    }
  }
  lambda.cov1 <- lambda.cov10 / sd.x 
  
  for (i in 1:dim.vill[1]){
    delta.village0[i] ~ dnorm(delta.village.cov10 %*% villcov[i, ] + delta.district0[district[i]], pow(sigma.district[district[i]],-2))
  }  
  delta.village <- (delta.village0 - mean(delta.village0[])) / sd.x  #Quantity of interest

  for (l in 1:dim.vill[2]) {
    delta.village.cov10[l] ~ dnorm(0, .01)
  }
  delta.village.cov1 <- delta.village.cov10 / sd.x
  
  for (i in 1:dim.dist[1]){
    delta.district0[i] ~ dnorm(delta.district.cov10 %*% distcov[i, ] + delta.province[province[i]], pow(sigma.province[province[i]],-2))
    sigma.district[i] ~ dunif(0, 10)
  }  
  delta.district <- (delta.district0 - mean(delta.district0[])) / sd.x  #Quantity of interest

  for (l in 1:dim.dist[2]) {
    delta.district.cov10[l] ~ dnorm(0, .01)
  }
  delta.district.cov1 <- delta.district.cov10 / sd.x

  for (i in 1:n.province){
    delta.province[i] ~ dnorm(0, .02)
    sigma.province[i] ~ dunif(0,10)
  }
  
  for (i in 1:dim.vill[1]){
    lambda.village0[i,1] <- 0
    for (k in 2:K){
      lambda.village0[i,k] ~ dnorm(theta.village.cov10[k, ] %*% villcov[i, ] + theta.district0[district[i],k], pow(psi.district[district[i],k], -2))
  }}  
  lambda.village <- lambda.village0 / sd.x #Quantity of interest

  for (l in 1:dim.vill[2]) {
    theta.village.cov10[1, l] <- 0
    for (k in 2:K) {
      theta.village.cov10[k, l] ~ dnorm(0, .01)
    }
  }
  
  theta.village.cov1 <- theta.village.cov10 / sd.x

  for (i in 1:dim.dist[1]){
    theta.district0[i,1] <- 0
    psi.district[i,1] <- 0
    for (k in 2:K){
      theta.district0[i,k] ~ dnorm(theta.district.cov10[k, ] %*% distcov[i, ] + theta.province0[province[i],k], pow(psi.province[province[i],k], -2))
      psi.district[i,k] ~ dunif(0, 10)
  }}  
  theta.district <- theta.district0 / sd.x #Quantity of interest

  
  for (l in 1:dim.dist[2]) {
    theta.district.cov10[1, l] <- 0
    for (k in 2:K) {
      theta.district.cov10[k, l] ~ dnorm(0, .01)
    }
  }
  
  theta.district.cov1 <- theta.district.cov10 / sd.x
  
  for (i in 1:n.province){
    theta.province0[i,1] <- 0 
    psi.province[i,1] <- 0
    for (k in 2:K){ 
      theta.province0[i,k] ~ dnorm(0, .02)
      psi.province[i,k] ~ dunif(0, 10)
  }}
  theta.province <- theta.province0 / sd.x #Quantity of interest

  omega[1] <- 0
  for (k in 2:K){
    omega[k] ~ dunif(0, 10)
  }

  sd.x <- sd(x)  #Essential to standardize results
  
  # For individual-level correlation between ideology and support for groups
  sd.s1 <- sd(s[,1])
  sd.s2 <- sd(s[,2])
  sd.s3 <- sd(s[,3])
  sd.s4 <- sd(s[,4])
  
  mean.x <- mean(x)
  mean.s1 <- mean(s[,1])
  mean.s2 <- mean(s[,2])
  mean.s3 <- mean(s[,3])
  mean.s4 <- mean(s[,4])

  cor.sx[1] <- sum((s[,1]-mean.s1)*(x[]-mean.x)) / ((N-1)*sd.s1*sd.x)
  cor.sx[2] <- sum((s[,2]-mean.s2)*(x[]-mean.x)) / ((N-1)*sd.s2*sd.x)
  cor.sx[3] <- sum((s[,3]-mean.s3)*(x[]-mean.x)) / ((N-1)*sd.s3*sd.x)
  cor.sx[4] <- sum((s[,4]-mean.s4)*(x[]-mean.x)) / ((N-1)*sd.s4*sd.x)

  ## Percent probability given to actual response (in-sample prediction rate)
  mean.correct <- mean(correct[,])
  
}
